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OS ■ Abstract 

, Proton transfer rates and mechanisms are studied in mesoscopic, 

Q>^ I liquid-state, molecular clusters. The proton transfer occurs in a proton- 

ion complex solvated by polar molecules comprising the cluster envi- 
ronment. The rates and mechanisms of the reaction are studied using 
both adiabatic and non-adiabatic molecular dynamics. For large molec- 
ular clusters, the proton-ion complex resides primarily on the surface 
of the cluster or one layer of solvent molecules inside the surface. The 
proton transfer occurs as the complex undergoes orientational fluctu- 
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' ations on the cluster surface or penetrates one solvent layer into the 

(— I ■ cluster leading to solvent configurations that favor the transfer. For 

Oh, smaller clusters the complex resides mostly on the surface of the clus- 

ter and proton transfer is observed only when the complex penetrates 
QJ [ the cluster and solvent configurations that favor the proton transfer 

are achieved. Quantitative information on the cluster reaction rate 
constants is also presented. 



X 

H , I. INTRODUCTION 

The rates and mechanisms of chemical reactions are influenced by the environ- 
ments in which they occur. Clusters with linear dimensions in the mesoscopic range 
are especially interesting reaction environments since their properties differ from 
either small molecular aggregates or bulk systems. Surface forces play an essential 
role in the dynamics as do molecular fluctuations. Since such clusters are likely 
involved in the reactive processes occuring in the atmosphere, the study of their 
reactive dynamics has practical as well as fundamental interest. 

A number of aspects of cluster reactions have been investigated recently. These 
investigations have examined the changes in the reactive dynamics that take place 
when a gas phase reaction is perturbed by clustering the reactive species with one 
or a few non-reactive or reactive molecules, reactive collisions where one or both 
of the reactive species are members of a cluster, cluster fragmentation reactions 
and ion- association reactions in water clusters. Proton transfer barriers in small 
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water clusters have been computed using density functional methods! and excited 
state proton transfer in large water clusters has been investigated experimentally .i 
The present article continues our earlier studyB of proton transfer in mesoscopic, 
polar molecular clusters. In contrast to this earlier investigation which focused 
on the activation free energy and its implications for the proton transfer reaction 
mechanism, here we consider the dynamics of the transfer process. We show that, 
within the mesoscopic domain, the mechanism is a function of the cluster size and 
describe the distinctive types of solvent dynamics that induce the proton transfer 
process. 

The model system is essentially the same as that in our free-energy study and 
similar to that in other studies of proton transfer in the bulk phase .i'i'Bi'i The 
system consists of a proton bound to a pair of A~ ions in a proton-ion complex, 
(AHA)~ . As in Refs.Di'i the distance between the two A~ ions in the complex is 
constrained to be 2.6A so that vibration of the A — A bond is not considered. The 
proton-ion intrinsic potential was constructed to model strongly hydrogen bonded 
systems.0 This intrinsic potential is given by: 

Vint = -I^IU' + \^2U' , (1) 

where = 1721.QK/A^, and ^2 = 12989.07^/1'^. The minima of the potential lie at 
±0.364A and the barrier height is about O.2/C5T. Here u is the proton coordinate 
along the ion internuclear axis measured relative to the center of the A — A bond 
in the ion pair. In the investigation of the dynamics of the proton transfer process 
the proton motion is restricted to the one-dimensional coordinate u. Earlier studies 
of proton transfer dynamics in the bulk phasei'i as well as our computations of the 
activation free energy for proton transfer in clusters^ allowed for the possibility of full 
three-dimensional motion of the proton. However, provided the intrinsic potential 
strongly confines the proton to the vicinity of the ion-pair internuclear axis, as was 
the case in these studies, the restriction to one-dimensional motion is an accurate 
representation of the dynamics. There are no fundamental difficulties in the removal 
of this restriction but the computation time increases. 

This complex is part of a cluster of diatomic molecules comprising the solvent. 
The diatomic solvent molecules are composed of two interaction sites with partial 
charges Za = ±0.52e. The bond length of the diatomic molecules was fixed at 2.0A 
giving a molecular dipole moment of yU = 5.0D. The interactions among solvent 
molecules, as well as those between the A" ions and the solvent, arise from site-site, 
6-12 Lennard- Jones (LJ) and Coulomb forces. The proton interacts with the solvent 
molecules via Coulomb forces. The parameters for the LJ potentials are described 
in.i 

The dynamics of the transfer process 

{A-H---A)- ^{A---H-A)- , (2) 

was investigated both by treating the proton degrees of freedom adiabatically and 
allowing for transitions among the protonic energy levels. Section |I| describes the 
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results obtained from the adiabatic simulations. This section begins with an outline 
of the simulation method and choice of reaction coordinate. The mechanism of 
the reaction that emerges from these calculations is then described, followed by 
a calculation of the rate constant of the reaction. Non-adiabatic dynamics is the 
subject of Sec. fT|. We first outline the surface hopping method!' which was used 
to carry out the non-adiabatic calculations. We then discuss the results and consider 
the modifications in the reaction rate and mechanism that arise as a result of non- 
adiabatic proton dynamics. Section |^ presents a summary and discussion of the 
results. 



II. ADIABATIC DYNAMICS 
A. Simulation method 



Let u be the proton coordinate and {R} the set of coordinates of all other 
classical particles in the system. The Schrodinger equation for the proton in the 
potential Vp depending on the fixed configuration {R} of the classical particles is 



^p(V„,m;{R})^„(«;{R}) 
e„({R})vI/„(w;{R}), 



+ Vp{u, {R}) 



2mp 



vl>„(«; {R}) = 



where Hp is the total proton Hamiltonian, mp is the mass of the proton and 27r^ 
is Planck's constant. The potential energy Vp is the sum of the intrinsic (|I|) and 
Coulomb Vc potentials, Vp = Vmt + The Coulomb potential Vc is given by 

i,a I -Tti.a " ttcM — Uttrel \ 

here a = 1, 2 labels the sites of any solvent molecule i, Rcm = (R/ + R'/7)/2 is the 
center of mass of the ion pair, tirei = (R-/ ~ ^ii) / I ~" R-// I is a unit vector 
directed along the ion-pair internuclear axis where R/ and R// are the positions of 
the two A~ ions. 

Classical particles with masses evolve according to Newton's equations of 
motion, 

mA = -Vr,K({R}) - Vr,(M/o({R}) I Hpim) I ^o({R})), (5) 

where the first term is force on Rj due to Coulomb and Lennard- Jones poten- 
tials that determine the solvent-solvent and solvent ion-pair interactions; the second 
term is the Hellmann-Feynman force that accounts for the action of the proton on 
the classical particles. Here | ^^oljR})) is the ket corresponding to the wavefunc- 
tion \E'o(w; {R}) in the u representation and Hp{{Il}) is the corresponding abstract 
Hamiltonian. 
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In order to solve @ the wave function \1/„(m; {R}) was expanded in a linear 
combination of localized Gaussian functions as 



^4u;{K}) = Y.c,{{K})Uu) , (6) 

with 



i=l 



Uu) = , (7) 

(2vra/) 

where fii and ai denote the position and width of the Gaussian functions. The 
proton-ion potential restricts the proton charge density to lie within the region be- 
tween the two ions and it remains localized along the internuclear axis joining the 
two ions. The intrinsic potential intersects the interionic axis at ±0.512A. Eighteen 
Gaussian functions were used to span the region between the A~ ions and the posi- 
tions of their maxima were located at equally-spaced points /Xj between —1.0 A and 
l.OA. The values of the widths af were taken to be 0.0225A for all basis functions. 
Use of the expansion given in results in a standard (non-orthogonal) eigenvalue 
problem : 

He = eSc , (8) 
where H is the Hamiltonian matrix with elements 

Hij = J (f)^{u)Hp(yu,u; {R})^j{u)du , (9) 

and S is the overlap matrix with elements 

Sij = J (f)i{u)(j)j{u)du . (10) 

The coefficients q satisfy the normalization condition 

n 

Y,CiSijCj = l. (11) 

The adiabatic dynamics calculation was performed as follows : The Schrodinger 
equation (|[) was solved for the ground state wave function and energy for a given 
configuration of classical particles. Note that the distance between the ions A~ 
is kept fixed so the time- dependent contribution to the potential arises from the 
positions of the solvent molecules in the cluster. Then, using the ground state wave 
function, the Hellman-Feynman forces given in (^ were computed and the classical 
equations of motion (|^) were integrated to yield a new classical configuration. The 
Verlet algorithmlll was used with time step of 5 x 10~^^s to integrate the classical 
equations. The constraints used to fix the intramolecular bond lengths were treated 
using the SHAKE algorithm.^ The constant temperature simulations were carried 
out using Nose dynamics.0 
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B. Rate constant 



One of the first steps in the study the rate of a reaction is the choice of a 
reaction coordinate. For adiabatic dynamics a natural choice for the proton transfer 
reaction coordinate is the expectation value of the position of the proton, Zp({R}) = 
(^o({R'}) I u I \[^o({R-}))- As we shall see, this reaction coordinate does provide 
a useful description of the proton transfer reaction dynamics. However, since it 
depends on the expansion coefficients of the ground state wave function which are 
known only numerically, it is convenient to seek an alternative reaction coordinate 
whose form is known analytically. Following earlierproton transfer studies! we use 
as a reaction coordinate the solvent polarizationE3'E3 given by, 

AE({R})^E..e(^-^). (12) 

Here s and s' are two points along the ion-pair axis in the reactant and product 
regions. We shall show that both the expectation value of the proton position and 
the solvent polarization provide equivalent descriptions of the reaction dynamics 
and either is a good reaction coordinate. 

Proton transfer in clusters is an activated process^ and if the free energy barrier 
is high enough a direct estimate of the reaction rate will require a long molecular 
dynamics trajectory. In this circumstance it is necessary estimate the reaction rate 
coefficient directly from the reactive-fiux correlation function.t3 Using the polariza- 
tion reaction coordinate, the (time-dependent) rate constant is given by 

_ {AE{0)6{AE{0) - AE^)e[AE{t) - AE^]) _ 

"^^^ ~ {e[AE{t) - AEt]) - ^ ■ ^^^^ 

where the angular brackets represent a canonical ensemble average, 6{AE{t)) is 
the Heaviside function and AE"^ is the value of the reaction coordinate at the bar- 
rier top {AE^ = for our symmetrical case). The rate constant is equal to the 
product of the transition state theory (TST) estimate of the rate constant and the 
transmission coefficient K{t). Using the Constrained- Reaction-Coordinate Dynamics 
(CRCD) ensembleEl the transmission coefficient is given by, 

^ ' {D-y^AEe{AE))^ ^ ' 

where (■ ■ ■)c is an ensemble average in the CRCD ensemble where the configura- 
tional distribution is taken from the AE = constrained dynamics but the velocity 
distribution is that for the system with no AE = constraint. The correction factor 
D that removes the bias generated by sampling initial configuration conditions from 
the constrained trajectories. For the polarization reaction coordinate D is simply 
given by, 
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2m ^ 



(15) 



where m = is the common mass of a solvent atom, i runs over the number of 
molecules and a = 1,2. The TST result may be computed from the expression 



TST 



(2vr/3)- 



-1/2 



{6{AE-AE^)) 



(16) 



C. Simulation results 

Calculations were performed for clusters of = 20 and N = 67 solvent molecules 
at temperatures of 200K and 260K, respectively. Under these conditions the clus- 
ters were in the liquid state. Evaporation did not occur on the time scale of the 
simulations, typically several nanoseconds. 

The fact that either the average proton position or the solvent polarization con- 
stitute acceptable reaction coordinates is demonstrated in Fig. |l| which shows that 
the time variations of both coordinates track the hops of the proton between the 
reactant and product configurations. Consequently we shall use both coordinates 
to provide insight into the reaction mechanism and for the computation of the rate 
constant. 



1. proton transfer mechanism 

The proton transfer mechanism consists of the description of the solvent dynam- 
ics in the course of the reaction (0). In our study of the proton transfer activation 
free energy we discussed the differences in the cluster structure when the proton was 
constrained to lie in the transition state or reactant (or product) regions.^ We may 
now describe the dynamical changes in the solvent that accomapany the transfer of 
the proton in the complex. 

We begin by confirming the picture of the cluster structure that emerged from 
the free energy study where a Feynman path intergral representation of the proton 
degrees of freedom was used and the centroidEl of the proton "polymer" was taken to 
be the reaction coordinate^ We consider the probability density, p±{z, r), for finding 
a positive or negative site on a solvent molecule at a point {z, r) in a cylindrical 
coordinate system centered on the A — A ion pair of the proton-ion complex with z 
directed along the A — A axis and r the radial coordinate in this cylindrical frame. 
The probability density is averaged over the angle variable. Rather than constructing 
this quantity by constraining the reaction coordinate to lie in the transition state or 
reactant (product) regions as was the case in Ref.,@ here we simply construct this 
quantity from a long unconstrained adiabatic molecular dynamics run and collect 
statistics for the probability density histogram only when the reaction coordinate 
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is found in the reactant or product configurations. Tlie results of this calculation 
are shown in Fig. ^ for a cluster with 20 solvent molecules and in Fig. ^ for a 67- 
molecule cluster. The insets in these figures show schematically the configuration of 
the proton-ion complex for the parts of the trajectory used to collect the statistics. 

The structural ordering in the cluster is evident in these figures. The cluster 
solvent molecules tend to strongly solvate the part of the proton-ion complex with 
the more exposed negative charge; i.e., the end of the complex which is less strongly 
bonded to the ion. This suggests that the complex tends to "float" on the 
surface of the cluster. When the ion is strongly bound to one A~ ion, this A — H 
dipole has a smaller dipole moment than that of a solvent molecule. Consequently, 
the solvent-solvent interactions are stronger than the interactions between a solvent 
molecule and this part of the proton-ion complex. These energetic arguments suggest 
that it may be favorable for this end of the proton-ion complex to reside on the 
surface of the cluster, a fact borne out by our simulation results. Of course, both 
entropic as well as energetic factors come into play in determining the structure of 
mixed clusters^ but here energetic factors seem to play a dominant role. 

There is also clear evidence of orientational order as can be seen from a compar- 
ison of the densities for the positive and negative sites. The same general picture 
applies for both the 20 and 67-molecules clusters. However, the 67-molecules clus- 
ter is large enough to support two solvent shells around the proton-ion complex 
and the presence of this second solvent shell has some important consequences for 
the dynamics which we shall describe below. The picture that emerges from these 
results is consistent with that from the earlier free energy studies: the position of 
the proton in the proton-ion complex has a strong influence on the structure of the 
cluster. 

Insight into the reaction mechanism can be obtained by examining the correla- 
tion between the time variation of the reaction coordinate and the solvent-complex 
dynamics. Let d, be the vector from the center of mass of the complex to that of 
the i^^ solvent molecule, = |(Ri,i + R-i,2) — R-cm = R-f^^ — R-cm, and d the 
vector from the center of mass of the complex to the center of mass of the solvent 
molecules in the cluster 

N N 

d = N-^Y. R-f - I^CM = ■ (17) 

i=l i=l 

Two quantities were used to gain insight into the nature of the solvent-complex 
dynamics: d =\ d \ and n* = n+ — n_ = YJi=i{(^{d ■ dj) — 6{—d ■ dj)) = — + 
2 X^ili ^(d ■ di). Here 6{x) is the Heaviside function and n+ and n_ are the number 
of solvent molecules with d ■ dj > and d ■ dj < 0, respectively. One expects small 
values of d when the complex is near the center of the cluster and large values when 
the complex is on the surface. Similarly, n* varies between n* = when the complex 
is symmetrically solvated in the center of the cluster and n* = N when it lies on the 
surface of a symmetrical cluster. It is possible to construct cluster configurations 
where the conditions are violated. Nevertheless, these quantities provide useful 
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indicators for solvent-complex dynamics and we have confirmed all aspects of our 
interpretations by direct visualization of the cluster dynamics. 

The results are shown in Figs. ^ and ^ In the case of a 20-molecule cluster 
one sees that the complex resides for long portions of time on the surface of the 
cluster (average distance d ~ 6) and long portions of time in the interior of the 
cluster (average distance c? ~ 2). If the complex is on the surface of the cluster 
then transitions rarely occur; however, if the complex makes an excursion into the 
interior (see the first large dip in n* and d in Fig. ^ then this excursion correlates 
with a proton transfer.Ey If the complex resides in the interior of the cluster, as is 
the case for times between 2.0 and 4.0 ns in the figure, then proton transfer events 
are frequent. 

The picture is somewhat different for the 67-molecule cluster. The proton-ion 
complex rarely penetrates deeply into the cluster. Some of the proton transfer events 
correlate with excursions of the complex into the cluster, but only one solvent layer 
deep and never far from the surface. However, even when the complex "floats" on 
the surface of the cluster there are frequent proton transfer events. Such events were 
not observed in the 20-molecule cluster case (although they may occur with lower 
frequency on longer time scales). The mechanism of these transfer events becomes 
clear from an examination of Fig. ^ which shows cluster configurations including the 
complex during one proton transfer event. Initially, the proton (black) is strongly 
bound to the "white" ion in the complex. The complex resides on the surface of the 
cluster with its "gray" end solvated in the cluster and it white end extending out of 
the cluster. In the course of time a fluctuation occurs that causes the complex to 
assume a configuration parallel to the surface so that there is a more nearly equal 
solvation of the two ends of the complex. This is a favorable configuration for proton 
transfer and transfer takes place around frame (d) of this figure. Once the proton 
transfer is complete, then the favorable complex configuration is for the now strongly 
hydrogen bonded gray end to protrude from the cluster and the weakly bound white 
end to lie within the cluster. 

Thus, the mechanism of the proton transfer depends on the size of the cluster. 
For smaller clusters fluctuations lead to penetration of the complex into the interior 
of the cluster where proton transfer is likely. For larger clusters transitions occur 
primarily by orientational motion of the complex on the surface of the cluster or when 
the complex makes shallow penetrations into the cluster. Deep penetrations of the 
complex into the cluster are rare. For the smaller 20-molecule cluster, presumably 
the orientational motion of the complex is restricted on the surface due to the 
larger surface forces. These features are borne out by an examination of the radial 
probability densities of the distance c?, Rc{d) shown in Fig. |^ (a) for a 20-molecule 
cluster and in (b) for a 67-molecule cluster. For the 20-molecule cluster there are 
two peaks, one corresponding to molecules in the interior of the cluster and the other 
for molecules on the surface. The results for the 67-molecule cluster show a broad 
single peak, with indications of some structure, corresponding to molecules on the 
surface and just within the cluster. 
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2. rate constant 

We now turn to a quantitative treatment of the reaction rate and compute the 
value of the rate constant. As described above, this calculation can be divided 
into two parts: the calculation of the TST estimate of the rate constant and the 
transmission coefficient. The product of these quantities yields the full rate constant. 

We begin with the computation of the TST rate constant using (|I^). The nu- 
merator of this expression is proportional to the probability density of finding the 
reaction coordinate at the barrier top {AE = 0) which in turn is related to the ex- 
ponential of the activation free energy. The free energy may be computed using the 
CRCD ensemble or estimated directly from a long unconstrained molecular dynam- 
ics trajectory if transitions are frequent. Figure § shows the free energy along the 
reaction coordinate determined in this way for a 67-molecule cluster along with the 
quadratic approximations to the results determined from fits around the potential 
minima. Note the parabolic nature of the free energy in the vicinities of the minima; 
in fact, the quadratic approximation is quite accurate over the entire range except 
at the barrier top where one expects the quadratic extrapolation to overestimate the 
barrier height. The direct simulation is also most subject to error at the barrier top 
where the probability density is smallest and here one expects the direct simulation 
to underestimate the barrier height. If one estimates the free energy barrier from the 
average of these two results one finds W{0) = 4.47kT for the 67-molecule cluster. 
Using (p!6[) one finds k^^^ = 0.04ps~^ = 1/25. Ops. The corresponding results for the 
20-molecule cluster are: W{0) = 7.75kT and = 0.025ps-i = l/40.0ps. 

The TST rate constant may also be estimated from the analytical formula 

iTST ^ ^0 -W(0)/kT 

27r 

where ujq is the frequency corresponding to the free energy minima. For comparison 
we may compute k^^^ using this formula along with the free energy determined 
from the proton position reaction coordinate instead of the solvent polarization 
coordinate. The results of this computation are: k'^^'^ = 0.036ps~^ = 1/28. Ops 
for the 67-molecule cluster and k'^^'^ = 0.014ps~^ = 1/69. 4ps for the 20-molecule 
cluster. The results for the 67-molecule cluster are in quite good agreement and 
the somewhat poorer results for the 20-molecule cluster arise from an inaccutate 
knowledge of the barrier height. 

The transmission coefficient was computed using (0). The following method 
was used to calculate the averages in the CRCD ensemble: Statistically independent 
classical configurations were selected every 10 ps from a long (1 ns in the case of 
20 solvent molecule cluster and 1.24 ns in the case of 67 solvent molecule cluster) 
constant temperature Nose molecular dynamics trajectory where the polarization 
reaction coordinate was constrained to lie at the transition state {AE = 0). Initial 
velocities were assigned according to the generalization of Boltzmann sampling for 
rigid diatomic molecules.0 For this ensemble of initial conditions, the constraint on 

9 



(18) 



the polarization coordinate was released and the trajectories were evolved forward in 
time using microcanonical molecular dynamics. This ensemble of trajectories could 
then be used to generate the averages needed to compute (0). 

Figures || (a) and (b), respectively, show the transmission coefficients as a func- 
tion of time for the 67 and 20-molecule clusters. /^From these graphs the transmission 
coefficients may be determined from the plateau values and one finds k = 0.4 and 
K = 0.5 for the 67 and 20-molecule clusters, respectively. The results in these figures 
show a rapid decay on a time scale which is less than a picosecond followed by a 
somewhat longer decay, of the order of a few picoseconds, to a plateau value. The 
time scale for the establishment of a plateau in K{t) is longer in the cluster environ- 
ment than in the bullJi for a similar but not identical intrinsic potential. This again 
signals different dynamics in the cluster compared to the bulk. 

The full rate constants determined from the product of k'^^^ and k are k = 
1/62. 5ps and k = 1/80. Ops for the 67 and 20-molecule clusters, respectively. Fi- 
nally, we may estimate the rate constant directly by simply monitoring the number 
of proton transfers in a long unconstrained adiabatic molecular dynamics trajectory. 
This procedure yields k = l/88ps for a 67- molecule cluster where the unconstrined 
trajectory had 77 proton transfer events. For a 20-molecule cluster where the tra- 
jectory had 29 transfer events the estimated rate is k = l/190ps. Once again the 
67-molecule results are in good agreement. However, for the smaller 20-molecule 
cluster, since the free energy barrier is higher and the quadratic approximation is 
valid over a smaller range it is more difficult to estimate the barrier height. In ad- 
dition, proton transfer is a more rare event so that the direct estimate of the rate is 
also subject to uncertainties. For these reasons the rate estimates are more variable 
for the 20-molecule cluster. 



III. NON-ADIABATIC DYNAMICS 
A. Simulation method 

The effect of non-adiabatic dynamics on the computation of the rate of the reac- 
tion was taken into account by using TuUy's surface-hopping, stochastic model0i 
which accounts for the possibility of quantum transitions in the dynamics of mixed 
quantum-classical systems. 

In this method a group of "classical trajectories" is considered. Each classical 
trajectory is evolved according to an equation of motion similar to (^) but with \E'o 
replaced by "^n, any of the adiabatic functions. The Hamiltonian (3) characterizes 
the quantum system. The wave function $(m,R, t) that describes the quantum 
mechanical state at time t, is expanded as a linear combination of the adiabatic 
states for the instantaneous "classical configuration" , 

{R}, t) = Y: C„W^n(w; {R}), (19) 

n 
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where C„ are complex-valued expansion coefficients. Substitution of the above equa- 
tion into Schrodinger equation yields the following equation for the evolution of the 
expansion coefficients: 

ihtu = Ck{Vkk - v^o) - ihY^c.-R ■ dky (20) 

j 

where 

Cj = Cjexpii fdt'Voo/^)- (21) 



Furthermore, using abstract notation 

\4,({R}) = (v&.({R}) I 4({R}) I vl>,({R})), (22) 
and the non-adiabatic coupling vector dfcj({R}) is defined as 

dfc,({R}) = i^kim) I V|K} I ^,({R})). (23) 



Equations ( |20D are integrated simultaneously with the classical equations of motion. 
Let denote by A the classical time step and 6 the quantum time step. At the end 
of each classical time step it is determined if a quantum transtion has taken place. 
According to the "fewest switches" algorithm the probability of switching from the 
current state k to all other states j during the time interval between t and t + A is 

_Mt + A)A 
- a,.(t + A)' 

where akj = CkCj and bjk = — 2_Re(a*^R ■ dkj). If Qkj is negative, it is set equal to 
zero. 

The simulation method is described in Ref.i Specifically, our simulations were 
carried out as folllows: Three quantum adiabatic states were used to describe the 
state of the proton. Initial conditions for a group of "classical trajectories" were 
determined. Statistically independent classical configurations were selected every 
10 ps from a long canonical run where the polarization reaction coordinate was con- 
strained at zero. Initial velocities were assigned according to the generalization of 
Boltzmann statistics for rigid diatomic molecules.ii Initially the total population was 
taken to be in the ground state, so Co = 1.0. For each of the "classical trajectories", 
initially two adiabatic steps were carried out to obtain the non-adiabatic coupling 
vector. Using the wave function the Hellmann-Feynman forces were computed. 
The classical equations (|) with \l/o replaced by were integrated using the RAT- 
TLE algorithms with a time step of A = 1 x 10~^ps. When the expectation value 
of proton position entered the transition region, —0.42 < {"^kit) \ u \ \I/fc(t)) < 0.42, 
starting from the previous time step, the integration step of the classical equations 
of motion was changed to A = 10^'^ ps and integration was continued for 150 time 
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steps. Equations (|20| ) were integrated using the Runge-Kutta method with time step 
S = 10~^ ps. At the end of each classical time step the switching probability was com- 
puted using (p^) to determine if a switch occured. If a switch occured conservation 
of the energy was satisfied by redistibution of the kinetic energy. Between transitions 
the coefficients Ck{t) evolve coherently. When (\l/fc(t- A) | u \ A)) < 0.591 < 

(^fc(t) I u I ^k{t)) or (^fc(t) I u I ^fe(t)) < -O.59A < {^k{t - A) I M I ^k{t - A)), 
so that the proton density enters the reactant or product regions, the coefficient for 
the current state was set equal to one, and all the other coefficients were set equal 
to zero. 

B. simulation results 

We may now examine the effects of transitions among the protonic states on the 
proton transfer dynamics. Figure |1^ shows a sample non-adiabatic trajectory for a 
67-molecule cluster. The lower panel in the figure is the proton position reaction 
coordinate computed using Zp = (\E'n|M|^'n), while the upper panel shows in which 
of the three states the proton lies. We note that the polarization coordinate is also a 
good reaction coordinate in the non-adiabatic case since it tracks the proton density 
changes. A number of features of these non-adiabatic trajectories are noteworthy. 

As expected, there is a strong correlation between transitions among the protonic 
energy states and the proton transfer events. The probability of a transition is 
large when the separation between the adiabatic energy levels is small and the 
separation is smallest in the transition region. Proton transfer does, of course, 
occur in the absence of transitions to the excited protonic states in accord with 
the predictions of adiabatic dynamics but there are substantial modifications to 
the simple adiabatic model because the separation between the ground and first 
adiabatic states is comparable to kT in the transition region. When the system is in 
an excited protonic state there is an increased probability of a proton transfer event. 
The proton density is more diffuse in the excited states than in the ground state. 
Also, when the system is in an excited protonic state and the reaction coordinate lies 
in the transition state region, the proton probability density is higher at extended 
spatial points where the proton ground state density is low. As a result, the proton 
density affects the solvent through Hellman-Feynman forces in ways that favour the 
transition state configurations. Consequently, when the proton is in the excited 
states re-crossings of the transition state region are more frequent than when it is in 
the ground state. This effect can be seen in the figure where the number of proton 
transfer events, and attempted proton transfer events, is larger when the system is 
in the excited states than when it is in the ground state. 

/^From these trajectories the picture of the proton transfer process is quite differ- 
ent from that for adiabatic dynamics since many proton transfer events are accom- 
panied by transitions to an excited protonic state. If the system makes a transition 
to an excited state it does not remain there long and quickly returns to the ground 
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state so that the majority of time is spent in the ground state configuration. ^Prom 
the simulation resuhs we estimated that the proton remains in excited states for 
only about 22 per cent of the total time. Thus, while in the course of passage from 
the reactant to the product regions the proton may make transitions into and out 
of an excited state, most of the dynamics is controlled by the ground state wave 
function. 

For non-adiabatic dynamics, the rate constant of the reaction was computed 
directly by dividing the time of the run by the number of proton transfers between 
product and reactant states. The simulation employed five trajectories (each 0.5 ns 
in duration) that followed microcanonical dynamics. The average temperature of 
the system was 260K. The rate constant was estimated to be A; = 0.017ps^^ (1/59.0 
ps) which is in good agreement with the adiabatic result. The fact that the rate in 
the non-adiabatic case is very close to that computed using adiabatic dynamics most 
likely arises from a cancellation of effects: when transitions to an excited protonic 
state occur and there is an increased probability of proton transfer, there is also 
an increased number of re- crossings in the transition region which lead to a small 
transmission coefficient and reduce the rate. 

IV. CONCLUSION 

Proton transfer rates and mechanisms differ significantly in the cluster and bulk 
environments. As in the bulk, solvation forces play an important role in determining 
the character of the cluster reaction but these sovation forces have a distinctive char- 
acter in the cluster. In the model investigated in this paper, the A — H end of the 
proton-ion complex, for instance in the configuration {A — H ■ ■ ■ A)~ ^ experiences in- 
teractions with the solvent dipolar molecules that are weaker than the solvent-solvent 
interactions. The • • • A~ end of the complex with the more exposed negative charge 
experiences strong charge-dipole interactions with the solvent molecules. These spe- 
cific forces are responsible for the tendency of the complex to reside on the surface 
of the cluster when the proton is in the reactant or product configurations and for 
the fact that the complex tends to be oriented normal to the surface with the • • • A~ 
in the cluster. 

There are strong fluctuations in these mesoscopic, hquid-state clusters. If the 
cluster is small enough, as in the 20-molecule clusters studied here, the complex 
resides predominantly on the surface of the cluster and reaction occurs only when 
it penetrates into the cluster. For larger clusters, e.g. the 67-molecule cluster, the 
complex again resides predominantly on the surface or one layer of solvent molecules 
inside the surface but now transfer is observed on the cluster surface as well. 

One may consider extensions of these ideas to other situations. Whenever there 
is strong charge separation in the course of reaction one might envisage scenarios 
like the one described above. However, the solvation forces could favor solvation of 
the complex in the interior of the cluster. This, in turn, would possibly give rise to 
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a transfer mechanism similar to that in the bulk phase. The interplay between the 
specific solvation forces, cluster size and the nature of the reactive species merits 
further investigation. 

The model for proton transfer considered here, that of a strongly hydrogen 
bonded system with negligible intrinsic barrier, favors the applicability of adiabatic 
dynamics. Nevertheless, we saw that although adiabatic dynamics does provide an 
estimate of the rate constant, transitions to excited states do occur and they have 
implications for the mechanism. Thus, it is interesting to investigate weakly hydro- 
gen bonded systems where non-adiabatic effects are even larger and could lead to a 
different picture of the proton transfer process. 
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FIGURES 



FIG. 1. Comparison of the proton position, Zp, and solvent polarization, AE (in units 
of lO^^^C/A), reaction coordinates. These reaction coordinates are plotted as a function 
of time for a = 20 molecule cluster. 

FIG. 2. Probability density of (a) the positive and (b) the negative solvent-ion sites 
in a cylindrical coordinate system centered on the proton-ion complex for a 20 molecule 
cluster. The z axis is along the A — A interionic axis and r is the distance to a solvent-ion 
site from the interionic axis of the A — A ion pair. The probability density was constructed 
by collecting configurations every 1.25 ps while the proton charge density is found between 
±0.7A and zbl.Ol. 

FIG. 3. Solvent ion probability densities. Same as Fig. ^ but for a 67-molecule cluster. 
The probability density was constructed by collecting configurations every 2.5 ps while 
the proton charge density is found between ±0.74 and itl.OA. 

FIG. 4. Trajectory of the proton position, Zp (top panel), the distance d between the 
center of mass of the cluster and the center of mass of the ion-pair (middle panel) and 
n* = n+ — n_ (bottom panel) as a function of time for a 20-molecule cluster. 

FIG. 5. Same as Fig. |^ but for a 67-molecule cluster. 

FIG. 6. Cluster configurations during a proton transfer event for a 67-molecule cluster. 
The configurations are separated by 2.5 ps. Time increases from panel (a) to panel (f). 

FIG. 7. Radial probability density Rc{d) = A-Kp{d)d^ versus d the distance between the 
center of mass of the ion pair and the center of mass of the solvent for (a) a 20-molecule 
cluster and (b) a 67-molecule cluster. The configurations used to construct these densities 
were collected every 1.25 ps for the 20-molecule cluster and every 2.5 ps for the 67-molecule 
cluster in the course of 4 ns unconstrained trajectories. 

FIG. 8. Free energy along the polarization coordinate for a 67-solvent molecule cluster. 
The open circles are the results from an unconstrained trajectory of 4 ns duration. The 
line is computed from the least-squares quadratic approximation to the numerical results 
at the free energy minima. 

FIG. 9. Transmission coefficient as a function of time for (a) a 67-solvent molecule 
cluster and (b) a 20-solvent molecule cluster. 
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FIG. 10. Non-adiabatic dynamics trajectory. The upper panel shows the protonic 
states and the lower panel gives the proton position reaction coordinate as a function of 
time. The error bars represent zt one standard deviation. 
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FIG. 11. Free energy along the polarization coordinate for a 67-solvent molecule clus- 
ter. The open circles are the results from an unconstrained trajectory of 4 ns duration. 
The line is computed from the least-squares quadratic approximation to the numerical 
results at the free energy minima. 
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FIG. 12. Radial probability density Rc{d) = A7rp{d)d'^ versus d the distance between 
the center of mass of the ion pair and the center of mass of the solvent for (a) a 20-molecule 
cluster and (b) a 67-molecule cluster. The configurations used to construct these densities 
were collected every 1.25 ps for the 20-molecule cluster and every 2.5 ps for the 67-molecule 
cluster in the course of 4 ns unconstrained trajectories. 



